source("/home/jc227089/evo-dispersal/KBGrad/KBGradfunctions.R")
setwd("/home/jc227089/SRE/KB/Adap/limits/data/")

args=(commandArgs(TRUE))

#evaluate the arguments
# input arguments will be b, h2D, fID
for(i in 1:length(args)) {
	 eval(parse(text=args[[i]]))
}

# modify the mother function to stop ztrace recording and break if population extent is too great
mother<-function(n, spX, K, lambda, b, ngens, Hmean, Dmean, h2H, h2D, VPH, VPD, evovar, stable.domain=spX){
  pop<-init.inds(n=n, spX=spX, Hmean=Hmean, Dmean=Dmean,
	h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPD)
	ztrace<-c()
  for(ii in 1:ngens){
    pop<-repro.disp(popmatrix=pop, K=K, lambda=lambda, b=b, 
    	h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPD, evovar=evovar, stable.domain=stable.domain)
    if (is.null(pop)) break
    if ((max(pop[,"X"])-min(pop[,"X"]))>300) break # stop looping if the population has begun to expand
    #if (ii>1000){
    #	ztrace<-rbind(ztrace, dzdx3(pop, K, lambda, b, stable.domain, ii))
    #}
  }
  out<-list(parameters=list(ngens=ii, K=K, lambda=lambda, b=b, n=n, spX=spX,
            ngens=ngens, Hmean=Hmean, Dmean=Dmean, h2H=h2H, h2D=h2D, VPH=VPH,
            VPD=VPD, evovar=evovar), pop=pop, ztrace=ztrace)
  out
}


stable.domain<-50
K<-150
lmda<-15
#b<-0.05
n<-2*stable.domain*K
spX<-stable.domain
ngens<-3000
Hmean<-0
Dmean<-log(4)
h2H<-0.5
#h2D<-h2H
VPH<-0.25
VPD<-0.05
nreps<-10

out<-vector(mode="list", length=nreps)
for (ii in 1:nreps){
	repl<-mother(n=n, spX=spX, K=K, lambda=lmda, b=b, ngens=ngens, Hmean=Hmean, 
		Dmean=Dmean, h2H=h2H, h2D=h2D, VPH=VPH, VPD=VPD, evovar=TRUE, 
		stable.domain=stable.domain)
	out[[ii]]<-repl
}

  	
save(out, file=paste("Hbyb_limits", fID, ".RData", sep=""))


	

